Statistical Mechanics of Interacting Run-and- Tumble Bacteria 
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We consider self-propelled particles undergoing run-and-tumble dynamics (as exhibited by E. coll) 
in one dimension. Building on previous analyses at drift-diffusion level for the one-particle density, 
we add both interactions and noise, enabling discussion of domain formation by 'self-trapping', and 
other collective phenomena. Mapping onto detailed-balance systems is possible in certain cases. 
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Several species of bacteria, including Escherichia coli, 
perform self-propulsion by a sequence of 'runs' - periods 
of almost straight-line motion at near-constant speed (v) 
- punctuated by sudden and rapid randomizations in di- 
rection, or 'tumbles', occurring stochastically with rate a. 
It is no surprise that the resulting class of random walk 
gives a diffusive relaxation of the number density at large 
scales [H- The resulting diffusion constant ~ v 2 /a, is 
vastly larger than that of non-swimming particles under- 
going pure thermal motion at room temperature. There- 
fore, apart from, e.g., the upper limit it imposes on the 
duration of a straight run (set by rotational diffusion), 
true Brownian motion can usually be ignored. 

Because bacterial diffusion is not thermal, the steady- 
state probability density cannot be written as p s oc 
exp[— H/ksT], with H a Hamiltonian, even for a single 
diffuser. The physicist's intuition can easily be led astray: 
for instance, Refs. 0,0,13] address models (of chemotaxis) 
comprising noninteracting particles in ID, with no exter- 
nal forces, but v(x),a(x) functions of position x. Instead 
of a uniform density, as would arise with any force-free 
detailed-balance dynamics, one finds p s (x) oc l/v(x) Q. 

Here we extend previous analyses of run-and-tumble 
motion to the many-particle level, addressing the roles 
of noise and interactions. These determine, for instance, 
the dynamic correlator of run-and-tumble bacteria, which 
is measurable by light scattering at low density [5( and 
at higher density, in principle, by particle-tracking mi- 
croscopy @. Additionally, particles for which v, a de- 
pend on the local density (either via thermodynamic in- 
teractions such as depletion or kinetic effects such as 
collision-induced tumbles) could show collective phenom- 
ena such as domain-formation or flocking. Such effects 
have previously been addressed within models where a 
self-propelled particle responds vectorially to the velocity 
of its neighbors, by direct sensing or passive hydrodynam- 
ics [1, 0, Below, we shall find, for run-and-tumble 
dynamics, similar effects in even simpler cases when only 
the speed of a particle is density-dependent. 

In making the transition from a single particle to many, 
most bacterial modelling approximates the number den- 
sity by a simple replacement p = Np 0, 0). But even 
for noninteracting particles, p (unlike p) is a fluctuating 



quantity, and a full statistical mechanics must compute 
noise terms for p. As seen below, these are not ad-hoc, 
but follow from the run-and-tumble dynamics directly. 

To allow relatively rigorous progress we work in 1-D 
throughout. For d > 1, although good descriptions exist 
at the one-particle level [3, we leave many-body ef- 
fects to future work. To avoid unwieldy equations for the 
probability flux J, we assume that, in units where run- 
and-tumble parameters such as v, a are 0(1) quantities, 
spatial gradients of these and of J are <C 1. (We show 
below that this was also implicit in [jj].) This restriction 
also judiciously avoids complicated memory effects aris- 
in g fro m the time-retarded response function of bacteria 
[Hlll|. that can lead to much reduced universality 12]. 

We start below with the microscopic dynamics of a 
single particle. We show how, on large scales, its prob- 
ability density is governed by a drift-diffusion equation, 
whose Langevin counterpart we then extend to describe 
the time evolution of the particle density in an assembly 
of interacting run-and-tumble particles. We then illus- 
trate some dramatic consequences of these interactions 
and finally address the role of external fields such as grav- 
ity. 

Let us define R(x,t),L(x,t) the probability of finding 
a single run-and-tumble particle at i,iina right- or left- 
moving state , respectively. With discrete run speeds 
vr,l and R <-> L, and interconversion rates cx.r,l/2 that 
depend on x, we have (with prime denoting d x ) 



R = -(v R R)' -a R R/2 + a L L/2 
L = (v L L)' + a R R/2-a L L/2 
P = -J' 



(1) 
(2) 
(3) 



Eq. ([3]) sums (1 1121) to relate the single-particle probability 
density p = R + L to its current J = v R R — vlL. Noise 
is not present in (|H2p . which already represent an exact 
master equation for the dynamics of one particle. 

The diffusive limit is found by a second time differenti- 
ation, elimination of R, L and their derivatives, and set- 
ting p — > 0] . The outcome can always be viewed as a 
time-local 'constitutive' relation J = J[p(x)] for use with 
the continuity equation ([3]). J, though linear in p, need 
not be local in space; in fact setting p = (which implies 
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J = 0) gives the differential equation (with 2v = vl ■ 
and 2a = qj, + olr): 



(1 
J 
V 

a 
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ti)J + &J' = J 

—Vp' + Vp; V = VRV^/a 
(ulvr - ctRVL)/2a - v(vrvl/v)' /a 
[vr(vr/v)' - v L (v L /v)']/2a 
(vr - v L )/a 



-Vr 

(4) 
(5) 
(G) 
(7) 
(8) 

An explicit nonlocal form J[p] is readily found from |4|. 
By neglecting J' in ([4]), we instead obtain the local form 
J = J /{I + £i) which is valid as long as V/V is itself 
small, a condition that must anyway hold if our diffusion- 
drift description is to avoid large p gradients in the flux- 
free steady state p s (x). The latter obeys V/V = (lnp s )', 
which follows exactly from (|H2p . Our local approxima- 
tion reduces to that of Schnitzer Q on further neglecting 
£i . While the latter maintains p s exactly, by retaining £i 
we extend this exactness to all steady states (J' = 0). 
We thus write as an optimal local approximation 



-J'\ 



J = -Dp 1 + Vp 



(9) 



with D = V/(l + £i) and V = V/(l + 

To analyse the more complex physics arising at many- 
body level, we first exactly recast ((9]) as an Ito-Langevin 
process for a trajectory Xi(t): 



Xi = A{xi) + C{x l )L i (t) 



(10) 



with Li(t) unit-variance Gaussian white noise, C 2 = 2D, 
A = V + dD/dxi, and C, A evaluated at the start of 



each time increment [14|. The latter (Ito) prescription 



eases the extension to many particles Xi(t) whose pa- 
rameters VR.L,ctR.L depend on position, not only ex- 
plicitly as in (fT0|) . but also implicitly via the particle 
density. The latter is formally constructed as p(x) = 
Xli Pi( x ) — J2i t>( x ~ x i)'i a coarse-grained version (com- 
patible with a gradient expansion) is found by choosing 
a smoothed function for pi(x). Within Ito's formulation 
[3], no further account need now be taken of the time- 
dependence of p: despite interactions, the random dis- 
placements C(xi)Li(t)St depend only on the preceding 
p(t) and are statistically independent of each other [lit . 

Accordingly, we can read Eq. (fl"0"|) as the Ito-Langevin 
equations for many interacting particles i = 1...N, with 
A = A(x, [p]) and C = C(x, [p]) [13]. From these, the Ito- 
Langevin equation for the collective density p{x) is found 
by standard procedures fl9| . as follows. Ito's theorem [l4[ 
states that for any function f(x{) of one variable 

f(xi) = (A + CL l )df/dx l +(C 2 /2)d 2 f/dx 2 (11) 
Piix^A + CL^f + Df']dx (12) 

Integrating p^|) by parts and using the identity f(xi) = 
J pi(x,t)f(x)dx gives 



Summing on i we obtain for the collective density [161 ] 
p = -(Ap)' + (Dp)"-J2(CL iPi )' = -J' c (14) 

i 

Jc = p(V + (5/Sp)'D) - Dp 1 + {2Dp) 1 / 2 K (15) 

with {A(x,t)A(x',t')) = 5(t-t')S(x-x') [II, [ll, Gil . In 
(fl"5|) . p is to be read as a coarse-grained, locally smooth 
field. Again by standard methods (but avoiding any ap- 
peal to detailed balance) [20] the Fokker-Planck equation 
for the many-body probability V[p] then follows as: 



V= dx 



5p{x) 



0, 



pV-Dd x p-Dp[d x 



Sp(x] 



V (16) 



Allowing that V, D are now functionals of p, we next 
seek conditions under which we can map (|16p , or equiva- 
lently (jTDJ) , onto a thermal system with detailed balance. 
(Clearly we require no macroscopic flux, (Jc) — 0.) In 
such a system, forces derive from an excess free energy 
T ex [p] ; diffusion and mobility matrices obey the Einstein 
relation Dij — pij; and steady state probability obeys 
V s [p] oc e"^W with T[p] = T ex + / p(\np - l)dx. (We 
set ksT — 1 without loss of generality. ) Since in (fT0|) the 
Li for different particles are independent, Dij is diagonal 
[l5j |. Hence the required condition is simply 

V([p],x)/D{[p\ 1 x) = -(6T ex [p}/6p(x))' (17) 

where the right hand side represents the force (i.e., excess 
chemical potential gradient) on a particle at x. 

For ID noninteracting particles, (JTTj) always holds, 
with Tex = J F(x)p(x)dx a,ndF(x) = J V(x')/D(x')dx' . 
For instance, when vl = vr = v(x) and qi = or ; = a(x), 
one has [3, 0] F(x) = lnv(x) so that the mean steady- 
state density obeys p s (x) — p s (0)v(0)/v(x). For interact- 
ing particles, existence of T ex is not generic even in ID, 
as the configuration space is A^ci-dimensional. Nonethe- 
less, some interesting cases do admit an T ex . In partic- 
ular, whenever V — (si(p))' and D = S2(p), with 
depending locally on p only, then V/ D satisfies (|17[) with 
Fex = / s 3 (p)dx and s 3 (p) obeying s 2 d 2 s 3 /dp 2 = dsi/dp. 

For example, consider a translationally invariant dy- 
namics in which vr^^olr^ at x depend on p(x), p'{x)..., 
with even parity: p(x) <-> p(—x) induces ccl(x), vl(x) <-> 
a R (-x),v R (-x). Then (5T ex {p]/5p{x))' = * + 
(ln(vRV L /v)y where * = (aRV L - a L VR) /2vrv L - On 
symmetry grounds, ^ must vanish when all odd deriva- 
tives of p do so; thus to leading order in gradients we can 
write \E' = (ip(p))' ■ Likewise vr^l — v oc ±p' so we can 
write {\xi(vrVl/v))' = 0- n v(p))' to leading order. Thus 
we recover T = J f(p)dx = J[p(lnp - 1) + f ex (p)]dx, 
with 



fex(p) = / [ip(u) +\nv(u)]du 



p i = -(Ap i )' + (Dp i )" -(LiCpi)' 



(13) 



(18) 



To leading order in gradients, this system is then equiv- 
alent to a system of Brownian particles at ksT — 1, with 
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mobility matrix Dij = SijV 2 (p(xi)) / a(p(xi)) , and a local 
free energy density /(p). Thus, according to mean-field 
theory [17] . whenever 



1000 



dp 2 



dp 



[ip(p) + \nv{p)} < 



(19) 



the system is locally unstable toward spinodal decompo- 
sition into domains of unequal density p. Also, it is glob- 
ally unstable to noise-induced (nucleated) phase separa- 
tion whenever p lies within a common-tangent construc- 
tion on f(p) [UllSl- In particular, for ijj = (left-right 
symmetry), any system with dv/dp < —v/p is liable, by 
the above reasoning, to undergo phase separation. Since 
one result of finite tumble duration is effectively to re- 
duce v [l3| , a strong enough tendency for the duration of 
tumbles to increase with density may have similar effects. 

In practice, of course, the very existence of a thermody- 
namic mapping in this (strictly ID) system ensures that 
the bulk phase separation predicted by mean-field theory 
is replaced by Poisson-distributed alternating domains of 
mean size cx e A . Here A is a domain- wall energy, fixed by 
gradient terms - if these remain compatible with the exis- 
tence of J-~ ex (which is not guaranteed) . To gain a first es- 
timate of such gradient terms (retaining ip = 0) we choose 
a model where, in a uniform system, v(p) — v Q e~ Xp with 
A a constant. Now we argue that the dependence v on 
p is somewhat nonlocal, sampling p on scales 7 of order 
the run distance, v/a, which is O(l) in our units. We 
therefore write for the nonuniform case v = v(p) where 
p = p + j 2 p" . (A linear term is forbidden by symmetry.) 
Thus SJ^ex/ '5p(x) = lnw = lnw D — X(p + 7 2 p") and 



f ex (p,p') =ph-iv + \[-p 2 + 1 2 p' 2 }/2 



(20) 



The resulting f(p) is very familiar [18[, and locally 
unstable for all p > 1/A. A stable dense phase is how- 
ever regained if v(p) saturates at large p. Suppose, e.g., 
v = v exp[— A<parctan(p/<p)], which falls as v Q e~ Xp be- 
fore approaching v sat — v D exp[— \(pir/2] at p <p. So 
long as ip > 2/ A, a window of phase separation is main- 
tained in mean field. Notably, some of the gradient terms 
arising from (the above form of) nonlocality now violate 
our thermodynamic mapping. However, if one ignores 
such violations, A is found to be large, and domain for- 
mation accordingly pronounced, whenever <p 3> 1. 

Although rigor is now exhausted, the physics seems 
reasonable: we know @, d, 0| that for an imposed 
v = v(x) particles accumulate in regions of low v. Thus 
with dv/dp sufficiently negative, 'self-trapping' of high- 
density, slow-moving domains can be expected. To in- 
vestigate whether this scenario arises, we have simulated 
Ecis. fTOf for the above conditions and indeed observed 
the predicted spinodal dynamics (FigfT]). Moreover, if 
we create a fully phase-separated initial state, this shows 
prolonged stability when the mean density lies between 
the predicted binodals; outside these, it collapses to uni- 
formity. Thus the self-trapping scenario appears valid 




t=0 



i = 100 



£=5000 



FIG. 1: Spinodal-like behavior of 6400 interacting particles 
within the region where d 2 f/dp 2 < 0. In the numerics, p 
is defined by convolution of p with a smooth function of fi- 
nite range ±1/2; we set v = 2.5, A = 0.01,ip = 250. The 
box length is 16, with periodic boundary conditions, and the 
mean density is 400. Dashed lines show the common-tangent 
densities (p = 60;800). 



despite violations of our thermodynamic mapping at gra- 
dient level. Eqs. (|14ll5p generalize obviously to d > 1, 
where they would lead to genuine phase separation un- 
der similar conditions. However, an adequate local ap- 
proximation relating A, C to run-and-tumble parameters 
remains to be established, and the range of validity of 
the diffusive limit is unknown. In d > 1 additional 
physics also enters, such as hydrodynamic interactions 
which are only partially accounted for by our use of a 
density-dependent velocity field [lol | . 

Finally, let us consider a translationally invariant sys- 
tem with no density dependence of a nor of v, but where 
vl.r are biased by some colloidal interaction H- m ^. That 
is, VR t L{%) = v ± vt, with vt = ^Pt^H^/ 8p{x))' and 
Pt the mobility (inverse friction) [7|. To first order in 
small Vt we then have V = vt and D — v 2 / a. The latter 
swamps any small thermal contribution (Dt = IzbTpt), 
so that (|10ll6p are, to this order, equivalent to a ther- 
modynamic system with H- lnt , but at enhanced temper- 
ature fcsTgff — D/ pt- Correlators such as S(q,t' — t) = 
(p q (t)p-q(t')) follow, although in many cases T e ff may 
be so large that these approach the noninteracting limit, 
S{q,t) = Nexp[-Dq 2 t] [21]. A similar expansion shows 
T c g- also to control sedimentation equilibrium under weak 
enough gravity (H = + rag J p{x)xdx). 

For the chosen p-independent v and a, this effec- 
tive temperature picture is intuitively clear and appeal- 
ing. Nonperturbatively, however, V = vt + (v T )'/a and 



V = (v 2 



'-)/a. Since vt — vt[p\, Eq. (fT7"|) no longer 



holds: the effective temperature concept breaks down as 
soon as the colloidal or gravitational interactions are non- 
infinitesimal. (This is true even within the gradient ex- 
pansion, which itself fails at vt/v ~ 1.) Perhaps instruc- 
tive is the exactly solved case of noninteracting sedimen- 
tation, where vt = —PTing- This gives in steady state 
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a Boltzmann-like exponential density, p s (x) = p(0)e~ KX , 
whose spatial decay rate k = —vtcc/(v 2 — Vy) is, how- 
ever, nonlinear in g [22| . On increasing g, the density 
profile collapses to zero height, not when g — > oo, but 
when vrig) — » v. Although cut off by other physics (such 
as true Brownian motion) this finite-<? singularity occurs 
because for \vt\ > v, L and R particles move in the 
same direction: a steady state of zero flux is clearly then 
impossible. Strong colloidal forces may likewise create 
absorbing states whose local 'escape velocity' exceeds v. 
This conclusion would alter if v represented not a fixed 
value, but the mean of a distribution extending (as in a 
truly thermal system) to unlimited speeds. Existence of 
a strict upper limit to the speed of self-propelled particles 
could thus be an important factor in their physics. 

In summary, by considering the passage from local run- 
and-tumble to drift-diffusion dynamics at large scales, we 
have elucidated the roles of both noise and interactions. 
When the mean run speed v is a sufficiently decreasing 
function of local density (or tumble-time t sufficiently 
increasing [HI), purely kinetic interactions could cause 
'self-trapping' of domains in ID, suggestive of bulk phase 
separation in higher dimensions. For particles interacting 
not by kinetic but by conventional thermodynamic forces 
(creating local drift velocities superposed on a density- 
independent run-and-tumble dynamics) such a mapping 
gives, perturbatively, an effective temperature set by the 
ratio of the run-and-tumble diffusivity to the thermody- 
namic mobility. However this mapping breaks down as 
soon as the drift velocities vt arising from the interac- 
tions become a significant fraction of the run-speed, v. 
For \vt\ > v, and in the absence of true Brownian mo- 
tion, absorbing states are possible. 

Because of the progress it allows, we have focussed 
above on those exceptional cases that, despite nonequi- 
librium interactions and noise, admit a thermodynamic 
mapping. Accordingly, large areas of parameter space 
remain unexplored; these could harbor many further in- 
teresting forms of collective nonequilibrium behavior. 

We thank R. Blythc, M. Evans, W. Poon and S. Ra- 
maswamy (discussions) and EPSRC EP/E030173 (fund- 
ing). MEC holds a Royal Society Research Professorship. 



[1] H. C. Berg, E. coli in Motion, Springer, NY 2004. 
[2] M. J. Schnitzer et al, Symp. Soc. Gen. Microbiol. 46, 
15-34 (1990). 

[3] M. J. Schnitzer, Phys. Rev. E 48, 2553-2568 (1993). 
[4] T. Hillen and H. G. Othmer, SIAM J. Appl. Math. 61, 

751-775 (2000); H. G. Othmer and T. Hillen, SIAM J. 

Appl. Math. 62, 1222-1250 (2002); R. Erban and H. G. 

Othmer, SIAM J. Appl. Math. 65, 361-391 (2005). 
[5] B. J. Berne and R. Pecora, Dynamic Light Scattering: 

with Applications to Chemistry, Biology and Physics, Sec. 

5.7, Wiley, New York (1975). 



[6] P. D. Frymier et al., Proc. Natl. Acad. Sci. USA 92, 6195- 
6199 (1995); G. V. Soni et al, Biophys. J. 84, 2634-2637 
(2003). 

[7] See, e.g., W. C. K. Poon, J. Phys. Cond. Mat., 14, R859- 
R880 (2002). 

[8] J. Toner, Y. H. Tu and S. Ramaswamy, Annals of Physics 

318, 170-224 (2005), and references therein. 
[9] A. Czirok, A.-L. Barabasi and T. Vicsek, Phys. Rev. Lett. 

82, 209-212 (1999); O. J. O'Loan and M. R. Evans, J. 

Phys. A. Math. Gen 32, L99-L105 (1999); G. Gregoire 

and H. Chate, Phys. Rev. Lett. 92, 025702 (2004). 
[10] A. Sokolov et al, Phys. Rev. Lett. 98, 158102 (2007); C. 

Dombrowski et al, Phys. Rev. Lett. 93, 098103 (2004); 

T. Ishikawa and T. J. Pedley, Phys. Rev. Lett. 100, 

088103 (2008). 

[11] P.-G. de Gennes, Eur. Biophys. J. 33, 691-693 (2004); 
D. A. Clark and L. C. Grant, Proc. Nat. Acad. Sci. USA 
102, 9150-9155 (2005). 

[12] Y. Kafri and R. A. da Silveira, Steady-state Chemotactic 
Response in E. coM larXiv:0710.4625t rl. eprint (2007). 

[13] We take tumbles as instantaneous, but a finite mean du- 
ration t is similar to setting (v,a) — > (v,a)/(l + ar). 
Indeed, for Poissonian tumble durations and vr,o.r, = 
vl,oll, the stated diffusive limit is thereby recovered, on 
neglect of all second-order time derivatives. 

[14] B. Oksendal, Stochastic Differential Equations, sixth ed., 
Springer, Berlin (2003). 

[15] Velocity-dependent forces, such as hydrodynamic inter- 
actions, would change this; see also [10(. 

[16] We have A(x„ [p]) = V(x„ [p]) + dD( Xi , [p])/dx„ whose 
second term is found by the chain rule as [D'(x) + 
(S/5p(x))' D(x)] x=Xi , and leads directly to the (5/Sp)'D 
term in (|15[) . This vanishes whenever D(xi, [p]) can be 
written D — D(xi,[p — pi(x)]), which excludes self- 
density interaction. 

[17] Omitting the noise term from (|15|l as in 0, [3[ amounts to 
a dynamic MFT (mean field theory), Jc[p] = Jc(x, [{p}], 
which loses the statistical mechanics of fluctuations and 
thus of activated phenomena such as nucleation. Static 
MFT instead finds the global extremum of T[p] , allowing 
access to nucleated phase separations [ilj ]. The thermo- 
dynamic mappings reported below include all noise terms 
and are therefore not contingent on either MFT. In some 
further cases, not explored here, such mappings can be 
made within dynamic MFT but are then violated by fluc- 
tuations. 

[18] P. M. Chaikin and T. C. Lubensky, Principles of Con- 
densed Matter Physics, Cambridge University Press, 
Cambridge (1995). 

[19] D. S. Dean, J. Phys. A Math. Gen., 29, L613-L617 
(1996); we extend this analysis here to the case where 
A, C are functionals of p. See also P. H. Chavanis 
[arXiv:0803.0263| 

[20] C. W. Gardiner, Handbook of Stochastic Methods, 3rd 
Edn., Springer- Verlag, Berlin (2004). See also, A. Lefevre 
and G. Biroli, J. Stat. Mech., P07024 (2007). 

[21] These results hold only at diffusive scales, qv/ot <C 1, 
whereas the opposite limit is more easily reached in scat- 
tering experiments Q. 

[22] Noninteracting particles in a quadratic potential, v T = 
— kx, obey p s (x) = pa\l — k 2 x 2 /v 2 \ a ^ 2k ~ 1 which is even 
further from Boltzmann form. 



